library(factoextra)
library(cluster)
diabetes <- read.csv("winequality-red.csv", header = TRUE, as.is = FALSE, sep = ';')
diabetes <- diabetes[, -1]
head(diabetes)
## class glucose insulin sspg
## 1 Normal 80 356 124
## 2 Normal 97 289 117
## 3 Normal 105 319 143
## 4 Normal 90 356 199
## 5 Normal 90 323 240
## 6 Normal 86 381 157
summary(diabetes)
## class glucose insulin sspg
## Chemical:36 Min. : 70 Min. : 45.0 Min. : 10.0
## Normal :76 1st Qu.: 90 1st Qu.: 352.0 1st Qu.:118.0
## Overt :33 Median : 97 Median : 403.0 Median :156.0
## Mean :122 Mean : 540.8 Mean :186.1
## 3rd Qu.:112 3rd Qu.: 558.0 3rd Qu.:221.0
## Max. :353 Max. :1568.0 Max. :748.0
Среднее и медиана отличается у всех признаков, ожидаем несимметричные распределения.
Построим matrix plot для того, чтобы увидеть особенности в наших данных.
library(lattice)
library(ggplot2)
library('GGally')
ggpairs(diabetes, title="correlogram", columns=c(2:4), upper = list(continuous = "points"), diag = list(continuous = "barDiag"))
Прологарифмируем признаки.
diabetes_l <- transform(diabetes, glucose=log(glucose), insulin=log(insulin))
names(diabetes_l)[names(diabetes_l) == 'glucose'] <- 'log_glucose'
names(diabetes_l)[names(diabetes_l) == 'insulin'] <- 'log_insulin'
ggpairs(diabetes_l, title="correlogram", columns=c(2:4), upper = list(continuous = "points"), diag = list(continuous = "barDiag"))
Удалим единичные outliers.
А также добавим раскраску по признаку “class”.
diabetes_lo <- diabetes_l
diabetes_lo[rownames(diabetes_lo)[diabetes_lo$log_insulin < 5],] <- NA
diabetes_lo <- na.omit(diabetes_lo)
ggpairs(diabetes_lo, title="correlogram", columns=c(2:4), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = class))
При дальнейшей кластеризации классифицирующий признак– “class” рассматриваться не будет.
library("FactoMineR")
dflo <- diabetes_lo[, -1]
res.pca <- PCA(dflo, scale.unit = TRUE, ncp = 2, graph= FALSE)
get_eigenvalue(res.pca)
## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 2.10928604 70.309535 70.30953
## Dim.2 0.83660334 27.886778 98.19631
## Dim.3 0.05411062 1.803687 100.00000
fviz_pca_biplot(res.pca, habillage=diabetes_lo$class)
Первые 2 ГК описывают 98% данных.
Рассмотрим 3 метода для определения числа кластеров:
Для точки \(i\) из кластера \(C_{I}\):
\(a(i)={\frac{1}{|C_{I}|-1}}\sum _{j\in C_{I},i\neq j}d(i,j)\) – среднее расстояние между i-той точкой и всеми остальными точками кластера.
\(b(i)=\min _{J\neq I}{\frac {1}{|C_{J}|}}\sum _{j\in C_{J}}d(i,j)\) – наименьшее среднее расстояние между i-той точкой и всеми точками любого другого кластера.
Теперь мы определяем силуэт (значение) одной точки данных:
\(s(i) = \begin{cases} 1-a(i)/b(i), & \mbox{if } a(i) < b(i) \\ 0, & \mbox{if } a(i) = b(i) \\ b(i)/a(i)-1, & \mbox{if } a(i) > b(i) \\ \end{cases}\)
Среднее s(i) по всем точками кластера - это мера того, насколько плотно сгруппированы все точки кластера. Таким образом, среднее s(i) по всем данным всего набора данных является мерой того, насколько правильно были кластеризованы данные.
Рассматривается характер изменения общего внутригруппового разброса с увеличением числа групп k. Объединив все n наблюдений в одну группу, мы имеем наибольшую внутрикластерную дисперсию, которая будет снижаться до 0 при \(k \rightarrow n\). На каком-то этапе снижение этой дисперсии замедляется - на графике это происходит в точке, называемой “локтем”.
Пусть \(E^∗_n{log(W^∗_k)}\) обозначает оценку средней дисперсии \(W^∗_k\), полученной бутстреп-методом, когда k кластеров образованы случайными наборами объектов из исходной выборки размером n.
Тогда статистика:
\(Gap_n(k)=E^∗_n{log(W^∗_k)}−log(W_k)\)
определяет отклонение наблюдаемой дисперсии \(W_k\) от ее ожидаемой величины при справедливости нулевой гипотезы о том, что исходные данные образуют только один кластер.
При сравнительном анализе последовательности значений \(Gap_n(k),k=2,…,K_{max}\) наибольшее значение статистики соответствует наиболее полезной группировке, дисперсия которой максимально меньше внутригрупповой дисперсии кластеров, собранных из случайных объектов исходной выборки.
fviz_nbclust(dflo, kmeans, method = "silhouette", k.max = 10) #"silhouette" (for average silhouette width), "wss" (for total within sum of square) and "gap_stat" (for gap statistics).
fviz_nbclust(dflo, kmeans, method = "wss", k.max = 10)
fviz_nbclust(dflo, kmeans, method = "gap_stat", k.max = 6)
Будем рассматривать деление на 2 и 3 класстерa.
set.seed(10)
km2 <- kmeans(dflo, centers = 2, nstart = 100)
km2$size
## [1] 117 27
km2$centers#/sapply(dflo, sd)
## log_glucose log_insulin sspg
## 1 4.74357 6.185327 139.9145
## 2 4.60564 6.145251 378.7037
km3 <- kmeans(dflo, centers = 3, nstart = 100)
km3$size
## [1] 52 9 83
km3$centers#/sapply(dflo, sd)
## log_glucose log_insulin sspg
## 1 4.562684 6.006067 239.4615
## 2 4.640043 6.213764 534.6667
## 3 4.823253 6.281513 112.4217
Кластеры получаются неравного размера (как и при рассмотрении классифицирующего признака, но там другое соотвношение индивидов в классах).
Деление на кластеры в плоскости первых 2 ГК и на ggpairs для 2 и 3 кластеров при использовании метода k-means:
fviz_cluster(km2, data = dflo)
fviz_cluster(km3, data = dflo)
ggpairs(dflo, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(km2$cluster)))
ggpairs(dflo, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(km3$cluster)))
Кластеризируем данные с меньшим числом признаков, для этого рассмотрим первые 2 ГК (на одну меньше, чем количество кластеров).
km3_pca <- kmeans(res.pca$ind$coord, centers = 3, nstart = 100)
km3_pca$size
## [1] 18 26 100
km3_pca$centers#/sapply(res.pca$ind$coord, sd)
## Dim.1 Dim.2
## 1 -0.8814062 1.8586132
## 2 2.8245783 0.1323704
## 3 -0.5757372 -0.3689667
fviz_cluster(km3_pca, data = res.pca$ind$coord)
#ggpairs(data.frame(res.pca$ind$coord), title="correlogram", columns=c(1:2), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(km3_pca$cluster)))
Деление также неравномерное.
Перейдем к следующему методу кластеризации– методу иерархической кластеризации с разными правилами объединения кластеров.
library(dplyr)
d <- dist(scale(dflo), method = "euclidean")
h_average <- hclust(d, method = "average")
plot(h_average, cex = 0.7, hang = -1) #агломеративная кластеризация, групповое среднее расстояние
cut_average <- cutree(h_average, k = 3)
df_cl <- mutate(dflo, cluster = cut_average)
ggpairs(df_cl, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(cluster)))
h_single <- hclust(d, method = "single")
plot(h_single, cex = 0.7, hang = -1) #агломеративная кластеризация, расстояние ближнего соседа (цепочки)
cut_single <- cutree(h_single, k = 3)
df_cl <- mutate(dflo, cluster = cut_single)
ggpairs(df_cl, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(cluster)))
Тут все совсем плохо ):
res.hc <- hclust(d, method = "complete") #агломеративная кластеризация, расстояние дальнего соседа (шарики)
plot(res.hc, cex = 0.7, hang = -1)
rect.hclust(res.hc, k = 3, border = 2:5)
cut_complete <- cutree(res.hc, k = 3)
df_cl <- mutate(dflo, cluster = cut_complete)
ggpairs(df_cl, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(cluster)))
Дивизионная кластеризация diana:
На шаге 0 все объекты объединены в один кластер.
На каждом шаге кластер делится, пока на шаге n - 1 все объекты данных не будут отделены друг от друга.
На каждом шаге делится кластер, назовем его R на два кластера A и B. Первоначально A равно R и B пуст. На первом этапе мы должны переместить один объект из A в B. Для каждого объекта i из A мы вычисляем среднее несходство ко всем другим объектам A:
\(d(i,A({i}))=\frac{1}{|A|-1}\sum_{i \neq j}d(i,j)\)
Объект i’, для которого уравнение выше достигает своего максимального значения, будет перемещено, поэтому положим
\(A_{new}=A_{old} \setminus \{i'\}\)
\(B_{new}=B_{old} \bigcup\{i'\}\)
На следующих этапах мы ищем другие точки для перехода из A в B. Пока A все еще содержит более одного объекта, мы вычисляем
\(d(i,A({i}))-d(i,B)=\frac{1}{|A|-1}\sum_{j \in A, i \neq j}d(i,j) - \frac{1}{|B|}\sum_{h \in B}d(i,h)\)
для каждого объекта i из A, и мы рассматриваем объект i’‘, который максимизирует эту величину. Когда максимальное значение вышеприведенного уравнения строго положительное, мы перемещаем i’’ от A до B, а затем просматриваем \(A_{new}\). С другой стороны, когда максимальное значение разности отрицательное или 0, мы останавливаем процесс, и разделение R на A и B завершено.
На каждом шаге делящего алгоритма мы также должны решить, какой кластер нужно разбить. Для этого вычислим диаметр
\(diam(Q)=max_{j \in Q, h \in Q} d(j,h)\)
для каждого кластера Q, доступного после предыдущего шага, и выбираем кластер, с наибольшим диаметром.
divisive.clust <- diana(as.matrix(d), diss = TRUE, keep.diss = TRUE) #дивизионная кластеризация
plot(divisive.clust, main = "Divisive", which.plots = 2, hang = -1)
cut_div <- cutree(divisive.clust, k = 3)
df_cl <- mutate(dflo, cluster = cut_div)
ggpairs(df_cl, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(cluster)))
Многие методы показывают похожие результаты. Если сравнивать с исходной классификацией, то методы кластеризации разбивают данные на другие группы.
divisive.clust <- diana(as.matrix(dist(scale(res.pca$ind$coord))), diss = TRUE, keep.diss = TRUE)
plot(divisive.clust, main = "Divisive", which.plots = 2, hang = -1)
cut_div <- cutree(divisive.clust, k = 3)
df_cl <- mutate(dflo, cluster = cut_div)
ggpairs(df_cl, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(cluster)))
Не слишком хорошо сработал метод понижения числа признаков.
library(mclust)
mc <- mclust::Mclust(dflo)
summary(mc) #первый символ относится к объему, второй - к форме, а третий - к ориентации. (E-- равный; V-- переменный; I-- расположение относительно осей координат)
## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 3
## components:
##
## log-likelihood n df BIC ICL
## -701.7836 144 29 -1547.692 -1559.786
##
## Clustering table:
## 1 2 3
## 33 83 28
head(mc$z) #каждому наблюдению назначается кластер с максимальной оцененной вероятностью.
## [,1] [,2] [,3]
## 1 4.352532e-03 0.9956475 1.496784e-14
## 2 1.344389e-07 0.9999999 8.739538e-39
## 3 8.770285e-07 0.9999991 4.940825e-32
## 4 2.443657e-03 0.9975563 2.713471e-15
## 5 5.391036e-04 0.9994609 2.013616e-19
## 6 9.513377e-03 0.9904866 7.771781e-12
par(mfrow = c(1, 2))
plot(mc, "classification")
plot(mc, "uncertainty")
fviz_cluster(mc, data = dflo)
Результаты похожи на исходную классификацию.
Максимальный BIC (среди всех моделей) у модели VVV (но у неё же и максимальное число параметров).
Попробуем выбрать другую модель:
plot(mc, "BIC") #Зависимость критерия BIC от числа кластеров для различных вариантов параметризации ковариационной матрицы
apply(mc$BIC, 2, which.max)
## EII VII EEI VEI EVI VVI EEE VEE EVE VVE EEV VEV EVV VVV
## 8 9 6 5 4 4 7 6 4 4 3 3 3 3
Кроме первых двух моделей BIC у остальных моделей близок, возьмем модель, кластеризующую данные на 3 кластера, но имеющюю меньше параметров.
mclustModelNames('EVV')
## $model
## [1] "EVV"
##
## $type
## [1] "ellipsoidal, equal volume"
mc1 <- mclust::Mclust(dflo, modelNames = 'EVV')
summary(mc1)
## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust EVV (ellipsoidal, equal volume) model with 3 components:
##
## log-likelihood n df BIC ICL
## -707.3937 144 27 -1548.972 -1560.778
##
## Clustering table:
## 1 2 3
## 25 88 31
head(mc1$z)
## [,1] [,2] [,3]
## 1 1.830318e-04 0.9998170 1.069044e-15
## 2 9.980166e-16 1.0000000 1.527052e-41
## 3 9.212325e-14 1.0000000 1.761841e-34
## 4 6.984441e-05 0.9999302 2.128655e-16
## 5 2.713809e-07 0.9999997 7.528674e-21
## 6 1.465955e-03 0.9985340 8.978124e-13
par(mfrow = c(1, 2))
plot(mc1, "classification")
plot(mc1, "uncertainty")
fviz_cluster(mc1, data = dflo)
Результаты немного отличаются, но все еще похожи на исходную классификацию.